function [x,output]=agentmaxrho_4agents_V2(xguess,index,ewsb,waut,bind,nobind);
global P beta n m S state agent gammagrid  Y scalefactor lowerbound rho table
objval = [];
output=[];
lb=[1.18 1.19 1.6];
ub=[1.19 1.20 1.81];
options = optimset('LargeScale','on','Display','iter','Algorithm','active-set');
[x,fval,exitflag]=fmincon(@objfun,xguess,[],[],[],[],lb,ub,@constrfun,options);
output(6)=exitflag;
function f = objfun(x);
i=index(1);
j=index(2);
consu(1:agent-1) = x(1:agent-1);
consu(agent)=Y(j)-sum(x(1:agent-1));
wint=interpolation_4agents(ewsb, index,x./Y(j));
for g=1:agent
output(g)=utility(rho,consu(g))+beta*wint(g);
end

objval=0;
for k=1:length(nobind)
objval = objval + (gammagrid(i,nobind(end))/gammagrid(i,nobind(k))-consu(nobind(end))/consu(nobind(k)))^2 ;
end
f=objval;
%
end
function [cineq,ceq]=constrfun(x);
i=index(1);
j=index(2);

consu(1:agent-1) = x(1:agent-1);
consu(agent)=Y(j)-sum(x(1:agent-1));

ceq=[];

for k=1:length(bind)
ceq = [ceq; ((utility(rho,S(j,bind(k)))+beta*P(j,:)*waut(:,bind(k))) - output(bind(k)))];
end
ceq
consu
cineq=[];
end
end